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ABSTRACT 

We modify the curvature condition for the existence of self-consistent scale-free discs, 
introduced by Zhao, Carollo & de Zeeuw. We survey the parameter space of the 
power-law discs, and show that the modified curvature condition is in harmony with 
the results of Schwarzschild's numerical orbit superposition method. We study the 
orbital structure of the power-law discs, and find a correlation between the popula- 
tion of centrophobic banana orbits and the non-self-consistency index. We apply the 
curvature condition to other families of scale-free elongated discs and find that it rules 
out a large range of power-law slopes and axis ratios. We generalize the condition, and 
apply it, to three-dimensional scale- free axisymmctric galaxy models. 

Key words: stellar dynamics - galaxies: kinematics and dynamics - galaxies: struc- 
ture - methods: analytical. 



1 INTRODUCTION 

A key problem in stellar dynamics is the determination of 
the range of axis ratios and density profiles for which triaxial 
galaxies can be in dynamical equilibrium (see, e.g., reviews 
by de Zeeuw 1996; Merritt 1999). Schwarzschild (1979, 1982) 
attacked this issue with his linear programming method, 
and constructed self-consistent triaxial galaxy models with 
constant density cores by numerical superposition of indi- 
vidual stellar orbits. Similar models with separable poten- 
tials were constructed by Statler (1987) and Hunter & de 
Zeeuw (1992). Gerhard & Binney (1985) showed that the or- 
bital structure in triaxial galaxies with central density cusps 
differs from those with cores. In particular, the box orbits 
are replaced by a host of minor orbit families and chaotic 
orbits (cf. Miralda-Escude & Schwarzschild 1989; Lees & 
Schwarzschild 1992). Based on this, Gerhard (1986) already 
suggested that cusped triaxial galaxies would evolve towards 
axisymmetry. Schwarzschild (1993) showed that moderately 
flattened scale-free triaxial models with an r~ 2 density law 
could be constructed numerically, but found that these mod- 
els evolved on long time scales. Similar results were obtained 
for non-scale- free models with steep central cusps (e.g., Mer- 
ritt & Fridman 1996; Merritt 1997; Siopis 1998). 

Although the general problem of the existence of self- 
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consistent triaxial, cuspy galaxies has not been solved, valu- 
able steps have been taken in the study of planar models that 
inherit many properties of the three-dimensional triaxial sys- 
tems. Kuijken (1993, hereafter K93) applied Schwarzschild's 
numerical construction method to elliptic discs with log- 
arithmic potentials, and found that self-consistent models 
are ruled out for axial ratios < 0.8. He traced the non- 
existence to the properties of the chaotic orbits, in line with 
the three-dimensional results of Schwarzschild (1993). Srid- 
har & Touma (1997a, hereafter STa) introduced a class of 
non-axisymmetric discs whose potentials are of Stackel form 
in parabolic coordinates. These discs admit an exact second 
integral of motion and host a continuous family of centro- 
phobic banana orbits. Syer & Zhao (1998, hereafter SZ98) 
showed by numerical means that none of the STa models 
can be self-consistent. 

Zhao, Carollo & de Zeeuw (1999, hereafter ZCZ) de- 
veloped a simple analytical tool for studying the self- 
consistency of scale-free discs. They compared the curva- 
ture of the model density with that of individual orbits near 
the major axis, and showed that this leads to a necessary 
(but not sufficient) condition for self-consistency which can 
be evaluated easily. They analysed a family of elliptic discs 
with power-law densities for a range of slopes, and showed 
that the curvature condition produces results in harmony 
with those of K93 and SZ98. We show here that while the 
ZCZ approach is sound, their condition for self-consistency 
contains an unfortunate error, which leads to incorrect re- 
sults for surface density distributions that do not fall off as 
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1/R. We present the correct condition for two-dimensional 
scale- free models, and compare the results with models con- 
structed by means of Schwarzschild's method. This not only 
validates the curvature condition, but also provides insight 
in the proper way to apply Schwarzschild's method. 

This paper is organized as follows. In we introduce 
scale-free elongated discs, describe how self-consistent discs 
can be constructed with Schwarzschild's method, and red- 
erive and correct the curvature condition of ZCZ. In §^ 
we study the parameter ranges of self-consistent power-law 
discs, and compare the curvature results with the numerical 
solutions. In §^ we apply the curvature condition to three 
other families of elongated discs, two of which are new, and 
confirm the results of K93 and SZ98. We extend the curva- 
ture condition to three-dimensional scale-free axisymmetric 
models in and apply it in §[| We summarize our conclu- 
sions in §[?]. We will study the triaxial case in a future paper. 



2 SCALE-FREE ELONGATED DISCS 

2.1 Surface-density, potential, and scaling 

We consider razor-thin power-law discs with a surface den- 
sity given by 



E(i?,0) = T R a ~ 1 S( ( / > ), 
and a corresponding gravitational potential 
f V R a P(<t>), a^O, 

\ V [2ln R + P(4>)\, a = 0, 



V(R,< 



(1) 



(2) 



where (7?, <f>) are the usual polar coordinates, and the func- 
tions S((f>) and P{(f>) are related by Poisson's equation. We 
are interested in models which are bisymmetric, and hence 
restrict ourselves to < <j> < it/2, where <j> — corresponds 
to the long axis. 

Since the potential and surface density functions of our 
discs are scale-free, the orbits at different energies are related 
by simple scaling factors in length and time (e.g., Richstone 
1980, hereafter R80). To find these factors, we use the equa- 
tions of motion: 



1 dV 

r ^ + 2R ^-rW 



(3) 



For a potential of the form (]2[) these equations are invariant 
under the scaling transformations 



R = kR, 



(4) 



where r and R are the scaled time and radius, respectively. 
This means that if a star with the coordinates [R(t) , <f)(t)] 
spends a time At in the angular sector Acj>, there exists 
another star with the coordinates [R(t), 4>{ t )] on the scaled 
orbit that spends the time 



(if 



At, 



(5) 



in the same angular sector. It follows that the velocity vector 
(vr, vj,) scales as (k^VR,k^v^>). 



2.2 Schwarzschild's method 

The numerical construction of the self-consistent discs de- 
fined by eqs (Q) and (^) is simplified by their scale-free 
nature. We follow the approach of R80 and K93, and at- 
tempt to reproduce the surface density on the unit circle, 
T(R = 1, <j>). We restrict ourselves to the first quadrant and 
divide the unit circle into N equally spaced azimuthal cells, 
with surface density E ce ii(i) = T(R = l,(j>i), (i = 1,..,N). 
We define E or b (i, j) as the mass contribution of the jth orbit 
(j — 1, .., M) to the ith cell. The model is self-consistent if 

M 

$^Eorb(M>0') =Eoeu(<), i = l,...,N, (6) 

3=1 

subject to the condition 



w(J) > 0, 



,Af. 



(7) 



The weights w(j) are a numerical representation of the 
phase-space distribution function that produces the surface 
density E in the potential V (cf. Schwarzschild 1979; K93). 

We take M 2> N, so that the system of equations (pi) 
and (Q) is overdetermined. A solution of (^) subject to the 
constraints (^) can be found by standard linear program- 
ming (LP) methods such as the simplex method (Press et 
al. 1992) or the non-negative least squares (NNLS) method 
(Pfenniger 1984). The NNLS method is useful when kine- 
matic constraints are included (e.g., Rix et al. 1997). We 
are interested only in reconstructing the surface density dis- 
tribution, and we employ the simplex method through min- 
imizing the non- self- consistency index Y defined in K93 as 



Y = 



JV 

NT ^ 



Ecoll(j) - 2J E olb (»,.?) 



(8) 



subject to the constraints w(j) > 0. The quantity Y must 
vanish for self-consistent discs. E is the mean surface density 
on the unit circle for < (f> < n/2. 

The calculation of E ce ii(i) is straightforward through 
eq. ([j]). However, care needs to be taken with the compu- 
tation of the orbital densities E or b(i,j). At some arbitrary 
time, stellar orbits do not necessarily cross the unit circle. 
But, some of their scaled, similar families do. The amount 
of mass that a star moving in the scaled orbit deposits in 
the azimuthal bin A<f> (on the unit circle) is proportional to 
Ar/A(f>. According to eq. (B), this is equal to R^~ x At/ A<p 
where we have set R = 1. Thus, having calculated the or- 
bital density at [R(t), <j>(t)], one can readily find the amount 
of mass contributed to the ith cell of unit circle through 
multiplying the obtained density by For a = 0, this 

factor is _R _1 as noted by R80 and K93. 

2.3 Curvature condition 

Consider an azimuthal element of length A<j> on the unit 
circle. The surface density per unit length of this element 
will be S{4>i)A(j> where 4>i is the angular position of the cen- 
troid of the ith element. The mass contribution of the jth 
orbit family to the assumed azimuthal cell is proportional to 



(At/T)jR 



a/2-1 



where (At/T)j is the fraction of the total 



integration time T that a test star of orbit family j spends 
in the angular sector {4>i, A(j>). The model is self-consistent 
when (cf. eq. |6|) 
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= 5(^)A0, i = l,.-., N, (9) 



M 

E 

3=1 

where M is the number of orbits, N is the number of az- 
imuthal cells on the unit circle, and the w(j) > are weights. 
We have taken time averages of orbital densities (denoted by 
angles) for all the times when an orbit returns to the same 
sector. We can rewrite (0) in the form 



M 



I (At/T)j R° 



/a-i 



A<f> S(<j>i) 



3=1 \ 



In the limit, as A<f> - 
also need to take M 



i N. 



(10) 



so that N ~ > oo and we generally 
> oo, one obtains 



If d 2 r/d0 2 of eq. ([TB]) does not have a zero, then the asso- 
ciated disc cannot be self-consistent. 

The condition for self-consistency derived by ZCZ ap- 
pears to be identical to eq. (|l5|), but their derivation did 
not take into account the scaling of the orbits, which caused 
them to use fi = /izcz = i?~ 7 S(0) and 7 = 7zcz = 1 — a. 
All the formulae given in ZCZ for the second derivative of 
F are still applicable — and are indeed identical to eqs ( |l5| ) 
to (j^|) — but they need to be evaluated with /j, and 7 given 
in eq. (Jl3[) . The definitions (and results!) agree only when 
a = 0, i.e., for discs with logarithmic potentials (i.e., fiat 
rotation curves) such as the models studied by K93. 



E 

3=1 



■w(j) 



R 



a/2- 



(11) 



This equation is the continuous version of ([]) in that the 
cells of configuration space are shrunk to zero size and in- 
stantaneous orbital densities are taken into account. The 
orbit families which occur in ( |ll"| ) can be regular or chaotic. 
Eq. ([n]) can be expressed in terms of the angular mo- 



mentum Jj = Rjyj 



]>>(j)<r) = i, r = 

3=1 

where we have defined 



1 



KRj 



\Jj\HRj 



7 = 1 + 



(12) 



(13) 



We assume S((j>) is smooth. Since Jj is also smooth by the 



equations of motion 
second derivative of ( 



we can follow ZCZ, and take the 
with respect to <j>j- This gives 



■■ 0, Vcfe € 



(14) 



3=1 

This relation is a necessary condition for F and will be sat- 
isfied if at least some orbit families have positive d 2 F/d<^ 2 
and others a negative value. If all d 2 T/d(f> 2 are either strictly 
positive or strictly negative, then the condition can only be 
satisfied by having positive as well as negative w(j), so that 
the model is not self-consistent. 

As pointed out by ZCZ, the evaluation of d 2 F/d<^ 2 
is simplified considerably near the major axis of the disc. 
Hence, we investigate (114) at cj)j — 0. We drop the subscript 
j of cj> and R but keep in mind that these variables are not 
the same for different orbits. Following ZCZ we obtain 

d 2 r 

d^ 2 
where 



^ = -£-[Qv-(l + 7-7A)], 



(15) 



and 



Kr = 



R 2 

ndV ' 

OR 



Ks = 



W) 2 

TfdV 
n dR 



The quantities Qv and Q M are defined as 



1 + 



n 8R 



<#>=o 



Q M = 1 + 



R 



dp 



<p=Q 



(16) 
(17) 

(18) 



3 THE POWER-LAW DISCS 

In order to compare the results of the curvature condition 
with numerically constructed discs, we now investigate a set 
of elongated discs with a very simple potential which makes 
it feasible to construct a large number of numerical models. 



3.1 Potential-density pair 

We take a potential of the form 
P(<j>) = {p 2 cos 2 <j> + sin 2 (f>) ~z , 



with 



(19) 



so that the equipotentials are similar concentric ellipses with 
axis ratio p. The angular dependence S((j>) of the associated 
surface density (|l|) is derived as a series expansion in Ap- 
pendix A, and given in eqs (Ai ) and ( A13). Figure la shows 
the surface density isocontours for p = 0.2 and a = 0.8. They 
are slightly dimpled near the minor axis (<j) = 7r/2), where 
S((f>) has a minimum. The dimples become stronger as p is 
decreased to zero. We have numerically evaluated S(ti/2) 
in this limit. Figure 16 shows that the result is positive for 
all a, and hence demonstrates that the surface density is 
positive for all values of < p < 1 and < a < 1. 



3.2 The curvature constraint 



The series expansion ( A£ , Ale ) for S(<p) can be used to show 
that the factor Q M in eq. (|16D is non-negative for < p < 1 
and < a < 1. As a result, < A m i n < A < A max < 00. 
Since 7 = l + a/2>0, we find that the second derivative of 
F will change sign when 



1 + 7 - 7A„ 



< 1 + 7 - 7A n 



(20) 



This is a necessary but not sufficient condition for self- 
consistency. This curvature condition imposes a constraint 
on the curvatures of V and fj, on the major axis. For 
Amin = 0, a non-trivial limit of (pOJ) is obtained as 



P > 



V2 + 



(21) 



which defines the boundary curve between the allowed and 
forbidden zones of the parameter space (solid curve in Fig- 
ure 2a). Specifically, p > 0.707 for a = and p > 0.633 
when a — 1. This corresponds to minimum axis ratios of 
the contours of constant surface density of 0.500 and 0.507, 
respectively. 
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Figure 1. Panel a shows the surface density isocontours of the power-law disc with p = 0.2 and a = 0.8. Panel b shows the lower bound 
of £ that occurs at <j> = tt/2 as p — > 0, as a function of a. 
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Figure 2. The (p, a) parameter space of power-law discs. For the marked points the reconstruction of the model density by Schwazrschild's 
numerical orbit superposition method was successful by assuming the convergence condition Y < ctjv, defined in equation ( ^3 ) (a) 
Squares correspond to N = M2 = 45 and filled circles correspond to N = M2 = 90 (b) Squares and filled circles correspond to the pairs 
(N = 45, M 2 = 180) and (TV = 140, M 2 = 560), respectively. In all cases we have set Mi = 80. The solid curve is the limit (0) of the 
region of non-existent models provided by the curvature condition Non-consistent models are to the left of the solid curve. The dotted 
line indicates the limit derived by application of the (erroneous) condition from ZCZ. The distribution of points in Panel b is discussed 
in §3.4 



3.3 Numerical models 

We now compare the results derived from the curva- 
ture condition with numerical models, constructed with 



Schwarzschild's method as outlined in §2.2. We divided the 
parameter space into 48 cells in the p-direction (from p — 0.3 
to p = 1) and 20 cells in the a-direction (from a = to 
a = 1). For each model we generated a library of M orbits 
containing Mi loops (flat tubes) and M2 orbits with zero 
initial velocities (along with their reflections with respect 
to the coordinate axes) dropped at M2 equally spaced az- 
imuths on the unit circle between (/> — and <f> = 7r/2 (these 



include fishes, bananas, pretzels, . . ., cf. Miralda-Escude & 
Schwarzschild 1989) . We compute the surface density of the 
ith cell at <j>i = — |) (1 < i < N), and uniformly 

distribute the initial positions of orbits (with zero initial 
velocities) on the unit circle, i.e., the jth orbit is dropped 
at <f>j = jn/[2(M 2 + 1)] where 1 < j < M 2 . In this way, 
initial positions of orbits will not coincide with the bound- 
aries of cells. Our loop orbit library consists of both thin 
and thick orbits, which are launched on the major axis with 
7/(0) = *(0) = 0, x(0) = 1 and y(0) > y min > 0. For 
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Figure 3. Variation of crjv versus a. The solid and dashed curves 
correspond to M2 = N = 90 and M2 = 4JV = 560, respec- 
tively. For each a, models having Y < <jn ar s accepted to be 
self-consistent. 



y(0) < y m in, the launched orbit will not belong to the 1:1 
resonant island. We find y m i n numerically. 

In order to calculate the orbital densities, we integrated 
the equations of motion for up to T = 120 galactic years us- 
ing the RK78 routine of Fehlberg (1968), which guarantees 
the preservation of energy with an accuracy of 10 -12 . For a 
regular orbit, T should ideally be the time by which the orbit 
becomes dense on its invariant torus, but this may require 
very long integrations for near-resonant and chaotic orbits. 
We therefore decided to follow Cretton et al. (1999), and to 
fix T for all orbits to a value that is sufficiently large for the 
remaining variations in the orbital density not to influence 
the main properties of the dynamical model. Experiments 
showed that for T > 120 the remaining 'noise' in the mod- 
eling procedure is dominated by the representation of phase 
space through a discrete grid, and not by the properties of 
the individual orbit densities. We therefore chose T = 120 
throughout. 

We ran our LP code for different choices of the num- 
ber N of angular cells, considered two values of the ratio 
M2/N, and took Mi = 80 loops in all cases. The results 
are illustrated in Figure 2, which shows all the points in 
the (p, a)-plane for which the numerical reconstruction of 
the model density was successful. For self-consistent discs 

Y should ideally become zero to within machine precision. 
However, our numerical computations show that for large 
values of p, the noise in the value of Y can have almost 
the same order- of- magnitude as Y itself. Thus, one needs 
a non-zero error threshold to pick up self-consistent mod- 
els. K93 defined a threshold value for non-self-consistency 
based on the isotropic distribution function of axisymmet- 
ric models, even though numerical experiments showed that 

Y is smaller for shallow axisymmetric cusps than for steep 
cusps. We therefore adopted the following scheme to find a 
credible threshold. For each value of a, we determine the 
function Y(p). We then filter Y(p) using the Savitzky-Golay 
smoothing filter (Press et al. 1992) to obtain the new distri- 
bution Y(p) from which one can calculate the noise distri- 
bution Ny(p) = Y(p) — Y(p). The models that we consider 
viable have Y(p) < <7jv where <tjv is the standard deviation 
of the noise distribution. It is defined by 

o N = - {Ny} 2 - (22) 



Figure 3 shows the variation of cr/v versus a for M2 — N = 
90 and M2 = 4iV = 560. According to this figure, shallow 
cusps (with (jm ~ 0.001) have more accurate Schwarzschild 
models than steep ones (with <tjv ~ 0.004). We note that 
most of our self-consistent models have Y(p) w 10 -8 , which 
is well within machine precision. For all successful models, 
noise fluctuations have been less than 1% of the global maxi- 
mum of Y . This indicates the acceptable quality of our mod- 
eling. 

Figure 2 a summarizes the results of our model com- 
putations for the case M2 = N, so that there is one zero- 
angular-momentum orbit dropped in each cell. We started 
our simulations with N = 45 (squares) and increased N un- 
til the set of self-consistent models (in the parameter space) 
converged (filled circles; N = 90). The lower bound of the 
allowed range of p increases as a tends to zero. Figure 26 
shows the same plane, but now for the case M2 = AN. It 
is evident that orbits with zero initial velocities play an im- 
portant role in the numerical reconstruction of the model 
density: by increasing M2 relative to N, the zone of admis- 
sible parameters grows considerably, a property that was 
also noted by K93. We discuss the implications of this nu- 
merical result, and the difference between panels a and b, in 
§3.4. 

In all models «90% of the loop orbits have zero weight, 
which indicates that they are less useful in the construction 
of self-consistent models. This is as expected (e.g., K93) be- 
cause the loop orbits of cuspy systems are universally anti- 
aligned with the potential and density isocontours, and they 
can hence only contribute part of the density (cf. de Zeeuw, 
Hunter & Schwarzschild 1987). As their individual orbital 
densities are smooth, a modest number of loops is sufficient 
to provide an accurate reconstruction of this part of the 
model density. Accordingly, the results do not depend sig- 
nificantly on the number Mi of loops in the orbit catalog, 
and our choice Mi = 80 is adequate for all cases. 

3.4 Implications for Schwarzschild's method 

Figure 2a shows that the possible region for self-consistent 
models predicted by the curvature condition agrees well 
with the results of the numerically constructed models when 
M2 = N: all of the numerically feasible models lie (well) in- 
side the possible zone when we use M2 = N. By contrast, 
when M2 = AN, many of the numerically acceptable solu- 
tions lie in the region forbidden by the curvature condition, 
even though the value of ctjv is smaller, so the condition on 
Y is more stringent. 

The choice M2 > TV is not consistent with the con- 
tinuum limit expressed by eq. (|ll|); as an azimuthal cell is 
shrunk to zero size, it becomes theoretically impossible to 
drop more than one orbit with zero initial velocity from that 
cell (it is still possible to launch an arbitrary number of or- 
bits with non-zero initial velocities from the same point). 
In the numerical simulations with M2 > N, more orbital 
'corners' are available per cell, which brings a degree of flex- 
ibility to the LP code that in turn results in more solutions 
that appear acceptable. Many of these are clearly ruled out 
by the analytic curvature condition. 

The requirement M2 = N also arises naturally in sepa- 
rable models, where the integral equation for the weights of 
the box orbits is solved analytically by assigning one box- 
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Figure 5. The variation of the relative population of banana 
orbits (Mb/2) versus p for several choices of a. The global max- 
imum of Mg/2 emerges in the neighborhood of p = 0.5 for all 
values of a. 



orbit corner to each cell (de Zeeuw, Hunter & Schwarzschild 
1987; Statler 1987). While box orbits have one corner in 
the first quadrant, regular boxlets may have more, so that 
the rule M2 = N may count some orbits more than once 
(K93). We conclude that in applications of Schwarzschild's 
method care should be taken that the number of orbits used 
is matched properly to the grid of cells on which the density 
is reproduced. 



3.5 Orbital structure 

To gain a better understanding of the origin of non-self- 
consistency, we constructed the Poincare maps of the orbits 
with zero initial velocities for several choices of a and p. We 
integrated the equations of motion in Cartesian coordinates, 
and sampled the phase space variables every time an orbit 
crosses the y-axis with x > 0. We have used the same library 
of orbits as in §3.3, and dropped M2 = 90 orbits from equally 
spaced azimuths between <f> = and (j> — tt/2. The energy 
E = \v 2 + V(R, 0) of all orbits was set to E = 1 and scaled 
orbits were used if needed, because orbits dropped from the 
unit circle do not all have the same energy. For each a, we 
generated the surfaces of section for two values of p: one 
in the forbidden zone and the other in the allowed zone 
of the parameter space explored by the curvature condition. 
The results are displayed in Figure 4. The non-self-consistent 
models (with p — 0.5) share an interesting feature: the bulk 
of phase space is occupied by banana orbits. By contrast, in 
models that are not ruled out by the curvature condition, 
the phase space is dominated by fishes and high-resonant 
orbits, while banana orbits are almost absent. This suggests 
that the non-self-consistency is related to the existence of 
centrophobic banana orbits that deposit much mass away 
from the major axis (as discussed also by, e.g., Pfenniger & 
de Zeeuw 1989; SZ98; ZCZ; Zhao 1999). 

To pursue this idea further, we computed the number 
Mb of banana orbits, using the information in the Poincare 



maps (in the y-y cross section; banana orbits dropped from 
the first quadrant never take y > 0). Figure 5 shows the ratio 
M B /2 = Mb /M2 (in percentage) versus p for several choices 
of a. M B /2 decreases when p — > 1 and has a maximum 
near p = 0.5. This result is not unexpected, for banana or- 
bits emerge when the frequencies of the x- and y-oscillations 
are commensurate in the ratio u x /ui y « 1/2. To first order, 
uj x /ujy equals the axial ratio p of the equipotential curves. 

Figure 6 shows the computed values of the non-self- 
consistency index Y (for N = M 2 = 90 and Mr = 80) 
versus p for the models with a — 0.8. The best-fitted poly- 
nomial curve has a similar variation as the corresponding 
Mb/2 curve in Figure 5. The models with 'banana-rich' or- 
bital structures (corresponding to p ~ 0.5) are maximally 
non-self-consistent. As p — » 1, Y decreases like M B /2 and 
converges to the required accuracy of self-consistent models. 

According to Figure 4, only very thin layers of chaotic 
orbits occur in the vicinity of hyperbolic fixed points (as also 
found by K93 for the logarithmic discs). Therefore, they do 
not seem to play an important role in the construction of the 
model density. This conclusion is supported by the following 
arguments. The orbits of a thin chaotic layer float around 
various resonant islands while they remain 'very close' to 
the outermost regular orbits (slow orbits) of these islands. 
It is possible to approximately generate such irregular orbits 
by a random superposition of their nearby slow orbits. Thus, 
in the limit when a chaotic orbit is fully mixed in the phase 
space, its orbital density approximately becomes equal to the 
density of the 'ensemble' of slow orbits. This means that a 
chaotic orbit enters Schwarzschild's method by adding (sub- 
tracting) a constant weight to (from) the weights of slow 
orbits. Our numerical computations confirm this point and 
reveal that our Schwarzschild models do not depend on the 
integration time of 'thin' chaotic orbits. The above reasoning 
is no longer valid for 'thick' chaotic layers, which, however, 
do not emerge in our models. We conclude that fishes and 
high-resonant orbits are the main building blocks of the self- 
consistent discs. 



4 OTHER ELONGATED DISCS 
4.1 Elliptic discs 

We now consider discs with elliptic surface density distribu- 
tions of the form nw, with 



S((p) = (p 2 cos 2 <j) + sin 2 



(23) 



The gravitational potential in the plane of the disc is of the 
form (§) with (e.g., Evans & de Zeeuw 1992 eq. 5.2) 



(p 2 cos 2 </>+sin 2 4>+u) 2 

[(l + u)(p 2 +u)]~^~ «5 



du, 



(24) 



where E and V are related by Vo = 27rGE p a /aB(i, ±=^) 
for a =£ 0. When a = we have V = 2GT, R F (l,p 2 , 0), 
with Rf the Carlson elliptic integral of the first kind (e.g., 
Press et al. 1992), and 



GE f ln(p 2 cos 2 0+sin 2 (j>+u) , 
P(<t>) = I ; du. 



Vo J [(1 + U )( p 2 +U ) M 1 





(25) 
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Figure 4. Poincare maps of the orbits with zero initial velocities in the power-law discs. Orbits are dropped from equally-spaced azimuths 
between <j> = and <j> = it/ 2. In all cases we have generated the surfaces of section for M2 = 90 orbits and have set the orbital energy 
E = \v 2 + V(R, <f>) equal to unity. 
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Figure 6. The variation of Y versus p for a = 0.8, Mi = 80 and 
N = M2 = 90. The solid line is the best fitted polynomial curve 
to the discrete numerical data shown by dots. Similar patterns 
occur for the other values of a. The non-self-consistency index Y 
has converged to zero with an accuracy better than crjv ~ 0.001 
(for p > 0.9 we arrived at Y « 10 — s ). 



The integrals ( pi} ) and (^BJ) for the 0-dependence of the po- 
tential need to be evaluated numerically. 

In order to apply the curvature condition, we evaluate 
Q M and Qv, and obtain 

_ 3 P 2 q + 2(1-q) ^ 



Q 



(2 + a)p 2 

„2x-P2(p,a) 



1 + (1-P 



Pi(p,a) 
where we have defined 



(26) 



P*(p,a) = 



'(p +^) 5_l (l + ^)"~dM, (i = l,2). (27) 



These two integrals can be evaluated in terms of beta and 
hypergeometric functions (Gradshteyn & Ryzhik 1980): 



P i (p,a)=p 2 - 2i B 



' 1 2i-l+q ' 
,2' 2 , 



l+ct 1 . 



i+f;l-p 2 ) . (28) 



They reduce to elliptic integrals when a — 0,1. 

A necessary condition for self-consistency is obtained 
from (|2o|) and (G^) by assuming A m i n = 0: 



d-p 2 )^H<i+f- 

Pi(p,a) 2 



(29) 



The boundary curve between the allowed and forbidden 
zones of ( |29| ) was computed numerically and is displayed 
in Figure 7 (solid line). The admissible axial ratios of the 
logarithmic discs have an infimum of p — 0.545. This result 
agrees with the findings of K93, whose numerical studies pre- 
dicted a lower bound of p ~ 0.8 for self-consistent models. 
As in the case of the power-law discs, the numerical solutions 
lie well inside the allowed region of parameter space. 



Figure 7. The parameter space of elliptic discs and projected 
power-law galaxies. The boundary curves between possible and 
non-consistent models are shown by solid and dashed lines for 
elliptic discs and projected power-law galaxies, respectively. Non- 
consistent models are to the left of boundary curves. The param- 
eter p is the axial ratio of isodensity contours of elliptic discs. 
For the projected power-law galaxies, however, the parameter p 
is related to the axial ratio of the surface density by eqs (B5), 
(B6) and (Bll). 



4.2 Projected power-law galaxies 

Projection of the density of scale-free triaxial models with 
potentials that are stratified on similar concentric ellipsoids 
(triaxial versions of the power-law galaxies of Evans 1994) 
provides elongated discs with surface-densities of the form 
(^) and associated potentials of the form (^). We have not 
found this potential-density pair in the literature, and give 
the ful l der ivation in Appendix B. The function S(c p) is given 
i n eq . (B1C), and the function P{4>) is given in eqs (BIS) and 
319,). We will refer to these discs as the projected power-law 



(gig). 

galaxies. 

On the major axis 
ies, we find 



= 0) of projected power-law galax- 
3a 2 p 4 + 6 - 6a + p 2 (9a-4-2a 2 ) 



= 1 + 



p 2 (2+a)(l+ap 2 ) 
1 — p 2 U2(p,a) 
p 2 Ui(p,a)' 



(30) 



where 
Ui(p,a) 



' 1 2i\-l+a s 



P 2 (l+a) 2 F 1 (Shi i;j+i + a ; l_pa) 



+ 



(2i-l) 2 Fi 



, 2 ; 1+1+ 2 ; 1 



V) 



(31) 



with i = 1, 2. One can easily verify that Q M > for all values 
of I a I < 1 and < p < 1. Therefore, the necessary condition 
for self-consistency becomes 



(l-p 2 )U 2 (p,a) 



<1+ 2- 



(32) 



p 2 Ui(p,a) 

The boundary curve between the allowed and forbidden 
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zones of inequality (g2|) is displayed in Figure 7 (dashed line) . 
In comparison with the power-law discs, the self-consistency 
of a larger fraction of the parameter space is ruled out for 
the projected power-law galaxies. 



4.3 Sridhar— Touma discs 

The separable discs introduced by STa have 

P(<f>) = [(1 + sin <t>) 1+a + (1 ~ sin <t>) 1+a ] , (33) 

where < a < 1. The associated density function S(<f>) is 
give n in eq. (4) of STa as a one-dimensional integral. From 
(p3J) we obtain Qy — 2 + a and Q M > 0. Upon substituting 
these results into (|2(j) we arrive at the requirement 

2+a<2+~ (34) 

which shows that self-consistent models are ruled out for all 
admissible values of a. This is identical to the conclusion 
drawn by SZ98, based on another method. The non-self- 
consistency of STa models can also be interpreted by our 
arguments of £3.5. All STa models, which have the fixed 
axial ratio p = 0.5 of equipotential curves (for all values of 
the cusp slope), are non-self-consistent because they host 
only a continuous family of centrophobic banana orbits. 



5 SCALE-FREE AXISYMMETRIC GALAXIES 

We now extend the curvature condition to scale-free axisym- 
metric models, in which the orbits can be described as two- 
dimensional motion in the meridional plane. 



5.1 Potential-density pairs 

Consider spherical polar coordinates (r, 8, <j>) where r is the 
radial distance from the origin, 8 is the colatitude and <j> is 
the azimuthal angle. We denote the momenta conjugate to 
(r,8,<f>) by (p r ,pe,P(p). In axisymmetric systems, the com- 
ponent Ptf, = L z of the angular momentum vector parallel to 
the symmetry axis is a conserved quantity. Furthermore, the 
variable cf> does not explicitly occur in the potential-density 
pairs, which take the following forms in scale- free models 



V 



V sgn(a)r a P(8), a ± 0, 
V [2 In r + P(8)], a = 0, 



p = por 



a<2. 



(35) 



5.2 Derivation of the curvature condition 

The curvature condition ( |l5| ) can be extended to three- 
dimensional axisymmetric models if we define F as 



I^IHoA')' 



where 



1/2 



(36) 



(37) 



is the modulus of the angular momentum and (cf. eq. [jyj) 



7 = 1+-. 



(38) 



The subscript j refers to the jth orbit family. The scaling 
properties of planar discs given in (Q) will be valid for three- 
dimensional models if we replace R with r. The smoothness 
condition for F becomes 



cFT 
d&j 



0. 



V6>, G 



o.- 



(39) 



where 6 = n/2 corresponds to the equatorial plane. This is 
identical to equation (|l4|), but now applies to axisymmetric 
systems. As before, we drop the subscript j and write 

1 d 2 r 1 d 2 L 



1 d 2 /i dl/ d/i 
Fd02 ~LlW~Jld(P + LdOjUe 



+2 



(40) 



We confine ourselves to the region near the equatorial plane. 
This simplifies the expressions considerably. For the models 
having reflection symmetry with respect to the equatorial 
plane, we have dp,/d6= dV/dd=0 at 8 = n/2. 

The first and second derivatives of L and p with respect 
to 8 can be obtained using the equations of motion 



r = p r , 



Pe 

= -ji Pe 



r =p r — — 
dV 



W + pI_ 

dr r 3 
Ll cos ( 



(41) 



(42) 



At 8 — 7r/2, eq. ( \i2\j implies pg = dpg/d8 ~ 0, which leads 
to r8 + 2r8 ~ near the equatorial plane. Furthermore, from 
@ we find that dL/d8 vanishes at 8 = tt/2. By using ^ 
we find 



1 d 2 L 
L d8 2 



1 8 2 V 

-2 02 d8 2 



(43) 



The derivatives of p: with respect to 8 can be expressed in 
terms of dr/d8 and d 2 r/d8 2 , which we evaluate by means 
of the equations of motion and the condition rd + 2r8 ~ 
which is valid near 8 = n/2. 

Combining all terms in expression (^), and carrying 
out some algebraic manipulations, then leads to 



d 2 r 



d8' 2 
where 

A = QftKe + (1 + l)K r + K<j,, 



v — (k + 7 — 



7A)]| 



with 

K r = 
and 



r 



{rOf 



1 L 2 Z 

or 



and the quantities Qv and are defined by 

8 2 V 



= 1 + 



r dV l e =2- 
dr 



1 + 



or 



(44) 



(45) 



(46) 



(47) 



(48) 
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These expressions resemble those for the elongated disc case, 
but contain the additional parameter re. 

A necessary condition for self-consistency follows from 
© as 

K + 7 — T^max < kQv < re + 7 - 7A min . (49) 

The greatest upper bound of (^) corresponds to A m i n = 0. 
Therefore, a non-trivial condition for self-consistency is 



Qv<l + ~- 

K 



6 APPLICATIONS 



(50) 



We now apply the curvature criterion to three families of 
scale-free axisymmetric galaxy models. 



6.1 The power-law galaxies 

The power-law galaxies of Evans (1993, 1994) have 

fa[m?(0)]^[2g 2 -m?(0)+ami(0)] , a # 0, 
S{6)=1 (51) 

{2[m\{e)Y 2 [2q 2 ~m 2 (9)], 



a = 0, 



and 

P(0) = 



[mm] 

lnm5(fl), 



a/2 



a = 0, 



(52) 



where Vo = iirGpo and we have defined m\(9) = q 2k sin 2 9 + 
cos 2 0. The density function p is positive for q 2 > (1 — a)/2. 
The parameter f3 introduced in Evans (1994) is equivalent 
to —a. 

From (^2|) we find Qv = 1/g 2 - Therefore, the condition 
( pp| ) reduces to 



, 2 > 



2k 



2k + 2 + a 



(53) 



This shows that the self-consistency of the power-law models 
is not ruled out near the equatorial plane. 

The above inequality has an interesting interpretation: 
since the quantity re is a dynamical property of orbit fam- 
ilies, as it involves not only the potential but also L z , we 
can investigate whether an orbit family is useful for the self- 
consistency of a given model. For example, k — 1 corre- 
sponds to orbits confined completely to a meridional plane 
(L z = 0) , and therefore, we can rule out self-consistent mod- 
els with 



v/2+T 



(54) 



using only the meridional orbits. T his i s the same result ob- 
tained for the power-law discs in 



3.2 



As re tends to zero, 
the amplitude of vertical motions decreases and more flat- 
tened models are supported. In the limit, the models with 
q — become self-consistent for re = 0. Thus, equatorial or- 
bits are essential building blocks of highly flattened models 
in general and limiting circular discs in particular. 



6.2 Scale-free spheroids 

Scale- free galaxy models with spheroidal density isocontours 
have widely been used in galactic dynamics (Qian et al. 1995, 
hereafter Q95). In this case 



S(e)= \ 
and 

1 (q 2 sin 2 9 + cos 2 9 + u) $ du 



P(9)={ 



■{q 2 + i 



(i+it) 

TrG Pn °f H<f si" 2 + cos 2 9 + u)du q = q 



(55) 



(56) 



f 

qV ° (l + u)(q2 + u)i 

where Vo = 2irGpoq a ~ 1 /a for a / and Vo = 
27tGpo-Rf(1, 1, q 2 )/q for a = 0. Rf is the Carlson elliptic 
integral of the first kind. 

Oblate models have < q < 1 and prolate models have 
q > 1. We restrict ourselves to the range — 1 < a < 1. 

The curvature parameter Qv (in the equatorial plane 
9 = 7r/2) is 



Qv = 1 + (1 - q )t—. r, 

Vi(q, a) 

where 

l/ i ( 9 ,Q)= 9 2 - 2l B(l,i + ^i) 2J Fi(l + f,l;i+ i ^;l-g 2 
for oblate models and 



V i {q,a) = q Mi B(l,i- 



(57) 

(58) 
(59) 



for prolate ones. By inserting expression ( ]57j ) in the condi- 
tion (po|), and noting that 7 = 1 + f-, we obtain the necessary 
condition for self-consistency as 



(l-q 2 )V 2 (q,a) 
Vx(q,a) 



K' + f) 



(60) 



Numerical computations show that the left-hand side of ( |60[ ) 
is positive for oblate models and negative for prolate ones. 
So, the condition (^) is fulfilled by the proper choices of 
< re < 1, and we cannot rule out any models. 

6.3 Oblate Sridhar— Touma models 

The oblate mass models described by Sridhar & Touma 
(1997b, hereafter STb) are the extensions of the flat STa 
models to three-dimensional space. These axisymmetric, 
cuspy models are integrable and admit an exact third inte- 
gral of motion. We have shown elsewhere (Jalali & de Zeeuw 
2001) that exact two- and three-integral distribution func- 
tions can be constructed for these models. It is useful to 
determine what the curvature condition gives. 
The STb models have 

S(9) = (a-cos(9)(l + cos<9) Q + (a + cos<9)(l-cos<9) Q , (61) 

and 



P(9) = (l+cos0) 1+Q! + (l-cos 



(62) 



where 1 < a < 2. It follows that Qv is given by Qv = 1 + a, 
which upon substitution in feOD yields the condition 



a < 



(63) 
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This is satisfied for a < 2 and all possible orbits correspond- 
ing to < k < 1, as expected based on the results of Jalali 
& de Zeeuw (2001). The phase space of STb models consists 
of a continuous family of reflected banana orbits all of which 
are useful for the construction of self-consistent equilibria. 



7 DISCUSSION 

Lynden-Bell (1962) showed that the distribution function 
of a stellar system depends on the phase-space coordinates 
through the isolating integrals of motion. However, such in- 
tegrals of motion are not available for most physical mod- 
els. In fact, there are few mathematical models of galaxies 
that have separable potentials. Finding the integrals of mo- 
tion is a formidable task for non-integrable systems. The 
existence of f(E,L z ) distribution functions for axisymmet- 
ric models (with E the orbital energy and L z the angular 
momentum parallel to the symmetry axis) has helped galac- 
tic dynamicists to develop more realistic galaxy models and 
extract their observable properties (e.g., Q95). But, in non- 
axisymmetric systems, only E is conserved, and we have to 
adopt other construction methods. Schwarzschild's (1979) 
paper was a major step in this way. His elegant numerical 
method allowed astronomers to attack the self-consistency 
problem of galaxy models from a different point of view, 
by constructing the coarse-grained distribution functions as 
weighted sums of the densities of individual stellar orbits 
without explicit knowledge of non-classical integrals of mo- 
tion. 

Schwarzschild's method does not require restrictions on 
the shape of the model or the functional form of the poten- 
tial, and can be applied to non-integrable models as well as 
integrable ones. The only drawback is the cost of the nu- 
merical computations, especially when the number of cells 
in the configuration space is large and the orbit library is 
rich. Although supercomputers can handle massive calcu- 
lations, a quick mathematical tool has always been attrac- 
tive. ZCZ developed one such tool, which is a simple mathe- 
matical test for self-consistency, and applied it to scale-free 
elongated discs. The idea derives from the fact that by in- 
creasing the number of orbits and cells in Schwarzschild's 
method, one can obtain a limiting analytical expression for 
the conditions of dynamical equilibrium. The smoothness 
condition for the quantity Y is such an interpretation of 
Schwarzschild's method. The curvature condition says that 
if the discrete numerical equations (subject to the positivity 
constraints on the weight functions) is solvable, a relation 
between the curvatures of the potential and density func- 
tions must be satisfied. This results in a necessary condition 
for equilibrium. It takes a particularly simple form near the 
major axis of the model where the surface density isocon- 
tours are strongly curved, and the reconstruction of the den- 
sity profile is most difficult, as it requires boxlets and chaotic 
orbits. We have corrected an error in the ZCZ derivation of 
the curvature condition, and also extended it to motion in 
axisymmetric potentials. 

The curvature condition is necessary but not sufficient, 
so we can use it to detect non-consistency, but cannot prove 
self-consistency. We can only state that in the allowed region 
of the parameter space, explored by the curvature condition, 
scale-free self-consistent models may be found. Our exper- 



iments on the power-law discs show that the results of the 
curvature condition (applied to the major axis of elongated 
discs) do indeed provide a boundary inside which the numer- 
ical solutions lie, and they furthermore have highlighted the 
need to be careful with the choice of the catalog of zero- 
angular momentum orbits when applying Schwarzschild's 
method. 

Our results reveal a correlation between the non-self- 
consistency index Y and the relative population of banana 
orbits Mb/2 ( as suspected in many of the earlier studies 
quoted). Maximally non-self-consistent power-law discs (for 
p « 0.5) have banana-rich orbital structures independent of 
the value of a. We speculate that this property is shared by 
all razor-thin discs whose potential functions allow for the 
existence of 1:2 resonant orbits, and conclude that it is of 
interest to search for more elongated discs without bananas 
(Jalali 1999; Sridhar & Touma 1999). 

In the case of three-dimensional axisymmetric models, 
we have employed our test near the equatorial plane where 
the circular orbits and equatorial rosettes overwhelm the ef- 
fect of other orbit families. The existence of the L z integral 
provides an extra quantity k that can be controlled to sat- 
isfy the necessary condition of self-consistency. This dynam- 
ical property of orbit families does not occur in the study 
of planar non-axisymmetric systems and is not expected to 
appear in triaxial systems. It is known that f(E,L z ) dis- 
tribution functions cannot be constructed for most prolate 
models (an upper bound exists for the allowed axial ratios 
of prolate models versus the cusp slope, e.g., Evans 1994, 
Q95), but the curvature condition does not reject the self- 
consistency of prolate galaxies. This is not a shortcoming, 
for the curvature condition is necessary but not sufficient. 

Our next goal is to extend the curvature condition to 
scale- free triaxial galaxies. The step to non-scale-free models 
remains as another challenging problem. 
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APPENDIX A: THE SURFACE DENSITY OF 
THE POWER-LAW DISCS 

We derive the surface density associated with the disc po- 
tential 



V{R, cf>) = V R a {p cos- 4 4> + sin 



< a < 1, (Al) 



with < p < 1. We introduce e = 1 — p and rewrite the disc 
potential in the form 



V= VoR a (l-2tcos<p + t 2 ) a/2 , 



(A2) 



where t = ecosci with |t| < 1 (we exclude the case e = 
1 from our calculations). The ei-dependence of V can now 
be expanded in terms of Gegenbauer polynomials as (see 
eq. [8.930] in Gradshteyn & Ryzhik 1980) 



V = VcTT J^C*^ (cos 0)t n , 
where A = —a/2 and 



(A3) 



Cr, (COS ( 



=£ 



r(A + fc)IYA + n - k) ,„, , , /A<N 



We now expand the term t n in a Fourier series, and write 



cos" tf> — S ' cos m<j>, 



where /3m = when n — m is odd and 



(A5) 



n — to ] , m/(l, 



(A6) 



J_ ( n 

I 2 "U/ 2 



m = 0, 



when n — m is even. Using expressions (A4) and (A5), the 
potential function becomes 



V = 



V R° 



oo n n 



(n) r(A + fc)r(A + n-fc) 
fe!(n - fc)![r(A)] 2 



n-0 m=0 k-0 

x [cos(2fc — n + m)(j) + cos(2fe — n — m)4>]. (A7) 



Now, we use eqs (Al) and (A5) of Syer & Tremaine (1996, 
hereafter ST96) to find the surface densities corresponding 
to single Fourier modes and superpose the results. Interest- 
ingly, the surface density becomes 



oo n n 



_ -Voir- 1 \-\- \- (n) r(A + fc)r(A + n-fc) 
2-n-a 2^2^2^ ePm fc!(n-fc)![r(A)] 2 



2ttG 

n=0 m=0 fe=0 
X [F(si) COS S\4> + F(S2) COS S2 



(A8) 



where 



F(s 



T( Sl /2 + a/2 + l)r(si/2 - a/2 + 1/2) 



r(s,/2 + a/2 + l/2)r( Sl /2-a/2) ' 
si = 2k — n + m, 

s 2 = 2k - n- m. (A9) 



The minus sign in front of ( A.8 ) is cancelled once the Gamma 
functions are evaluated. 

The case a = 0, which corresponds to logarithmic po- 
tentials, is special and needs a different treatment. In this 
case we have 



V(R,cf>) = V [2m_R + ln(l - 2£cos0 + t 2 )] 



(A10) 



The second term in ( |A10| ) can be expanded as the series 
(Gradshteyn & Ryzhik 1980, eq. [1.514]) 



ln(l - 2t cos 4> + t 2 ) = -2 V 

^ — ' n 



(All) 



which after the substitution t = e cos <f) and using ln(l — z) = 
— X^^Li z " l n ' S ives the potential as 



V 



OO 

Vo { 2 ln i 1 - 1) + 2 ln R - £ 2^ cos 2n4> 



££ 



s (n + m) <(> + cos (n-m)4>] | . (A12) 
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By dropping the constant term in V and using eq. (A5) of 
ST96, we obtain the surface density as 



oo 



with k — 2n, z = n + m and j = n — m. 
A simple representation for E is 

1 °° V 

£(i?, 0) = — (a + ^^ae cos AM, a = 



(A13) 



(A14) 



where the Fourier coefficients at [t > 1) are given by 



= 
with 



7rG ^ 



n r 



r + E 



2n 



1, r = s, 
0, r/s. 



(A15) 



(A16) 



From the definition of ftm , we conclude that ai — when 
£ is odd. Using the rule F(z + 1) = zT(z) and eq. ( A5) o f 
ST96, the potential function can be regenerated from ( A14 ) 
as 



V = 2-kGcio \nR- 2nG — ' cos £4>, 



(A17) 



which is identical to the result obtained by K93 for aligned 
discs (the coefficients at of K93 are computed for an elliptic 
disc with flat rotation curve). 



APPENDIX B: THE PROJECTED 
POWER-LAW GALAXIES 

We compute the projected surface density of the scale-free 
triaxial power-law models, for arbitrary direction of viewing, 
and show that the resulting surface density has a simple 
form. We consider it as the surface density of an elongated 
scale-free disc, and show that the gravitational potential in 
this disc can be evaluated explicitly. This then defines a new 
family of elongated discs which we call the projected power- 
law galaxies. 

Bl Scale-free triaxial power-law models 

Consider scale-free triaxial power-law models with gravita- 
tional potential 

/ 2 2 \ 

!„»ln^» + a = , 

V(x,y,z)={ ' (Bl) 

The associated density distribution is (Binney 1981; de 
Zeeuw & Pfenniger 1988) 

- J Ax 2 + By 2 + Cz 2 



47rG (x 2 +p~ 2 y 2 + q~ 2 z 



(B2) 



with 

A = 

B = 



n2 „2 



1 - 



(B3) 



This is the triaxial generalization of the axisymmetric power- 
law galaxies of Evans (1993, 1994). The surfaces of constant 
density are not ellipsoidal and become dimpled near the 
short axis when q is small (e.g., Evans, Carollo & de Zeeuw 
2000). The density distribution ( [B2]) is non-negative as long 
as q 2 +p 2 q 2 — (1 — a)p 2 > 0. This is a non-trivial constraint 
except in the limit a — > 1. 



B2 Projected surface density 



Now observe the density (B2) from an arbitrary direction 



defined by the standard spherical polar coordinates (0,<j>). 
Define a coordinate system (x" , y" , z") such that the z" axis 
lies along the line-of-sight and the x" axis lies in the (x,y)- 
plane. The relations between (x,y,z) and (x" ,y" , z") are 
given in, e.g., de Zeeuw & Franx (1989). Direct integration 
shows that E(x" ,y") is an elementary function, given by 



E(x" 



where 

• 2 

c\ = sin 
c 2 = 2(1 

2 

C3 = COS 

and 



pqB 



T i 

,2 > 2 ' 



4ttG 
a\x 



(cic 3 

//2 



a2x"y" + a 3 y" 2 



(cix 



111 



c 2 x"y" + c- A y 



» 2 1 



2 2 
i + p COS ( 



P 

b cos 



sm <p cos (p cos ( 

2 n i 2.2, 

p sm <pi 



(B4) 



(B5) 



+ q 2 sin 2 ( 



ffli = C1C3 
a 2 = QC2(ci+c 3 ) 
a3 = C1C3 



12, / 2 , 1 2\ 



(B6) 



1 2 / 2 , 1 2\ 

jc 2 + a c 3 + z c 2 



The surface density is scale-free, and has biaxial symmetry. 

Now transform to polar coordinates (R, $") defined by 
(x",y") = (i?cos$",i?sin$"). Then the major axis of E 
lies at position angle $" = $*, with 



tan2$* = 



c 2 



0,2 



(B7) 

C3 — ci 03 — ai 

This is the expected result: as the potential is stratified on 
ellipsoids, its projection is stratified on similar concentric el- 
lipses at the fixed position angle $* (Stark 1977), and this is 
then also the position angle of the projected surface density. 

When we rotate the coordinate system to (R, $) where 
$ = <!>" — $*, we obtain 

E(R,$) = X R a - 1 S(ji), (B8) 
where 



En = 



)pqP 



4ttG 



(Cl +C 3 - A) 2+ 2 



(B9) 



© 0000 RAS, MNRAS 000, 000-000 



14 M. A. Jalali and P. T. de Zeeuw 



and 



(l-a)P 2 +a(l + P 2 )fi 



with 

r,2 Ci + C 3 



V(ci - c 3 ) 2 



(BIO) 



(Bll) 



ci + c 3 + A ' 

and we have defined the abbreviation p = P 2 cos 2 $ + sin 2 $ 
which will be useful below. 

Unlike p, the projected surface density E is everywhere 
positive. The angular dependence of the projected surface 
density has the same form for all directions of observation, 
but the value of P is a function of the intrinsic axis ratios 
p and q and the viewing angles 6 and 4>. The isophotes are 
slightly oval. An example is shown in Figure 4 of Evans, 
Carollo & de Zeeuw (2000) . The axis ratio P' of the surface 
density, defined as S(x' , 0) = S(0, P'x') is given by 

(B12) 



P' = 



[{1 + a)(l + ap 



2\]l/(l- a ) 



which reduces to P when a = 0. When the model is observed 
down the z-axis, we have <j> is arbitrary, 9 = 0, and P — p. 



when a = 0. Here Rf is the Carlson elliptic integral of the 
first kind (e.g., Press et al. 1992), and we have corrected 
some typographical errors in expressions (5.10) and (5.11) 
of Evans & de Zeeuw (1992). 



For the speci fic ca se with S(fj,) given in eq. (B10), the 
integral equation (B13) can be solved to give (with P = p) 

w(t) = a(l+p 2 )8(T) + 2 P 2 8'{t), (B17) 

where <5(r) is the Dirac delta function and S'(t) is its deriva- 
tive (see, e.g., Morse & Feshbach 1953, p. 837 or Andrews 
& Shivam oggi 1986, p. 35 ). The delta functions reduce the 



integrals ( B16 ) and ( B14 ) to single integrations. We find 

oo 

2irGT, p 2+o 'R a f {u+p)%du 



V(R,fl) 



aB (i!-f) Juh[(i+u){p 2 +u)]h+% 

l+a 



X 



p 2 +u 



+ 



o 

l+a 



- u 



U + fl 



(B18) 



for a/fl. For a = the coefficients A and B in eq. ( Bl. 
are given by 

A = ^p 2 GZ {R D (0, l,p 2 )+R D (0,p 2 , 1)], 



B3 Potential in the disc 



Now con sider surface densities of the form ( B4 ) with S(fi) 
given by (BIO) where without loss of generality we take P = 
p and <& = <j) (see above). We take Eo given, but note that 
( |B9[ ) could be used to relate it to a three-dimensional model. 

The potential V(x, y) in the plane of the disc can be 
evaluated with the formalism developed by Evans &i de 
Zeeuw (1992, §5). They decompose a given biaxial scale-free 
disc into a weighted integral of constituent discs with ellip- 
tic surface densities of different axis ratios. As the potential 
of the constituent discs is known from classical theory, their 
weighted integral provides the potential of our given disc. 
Specifically, if we can find the function w(t) such that 



Bifi) = 2p 2 GE 



\n(u + n)du 



u+p 2 u+lj ^/u(u+p 2 )(u+l) 



-^P 2 GE 0j Rj(0,p 2 ,1,m), 



(B19) 



where Rd and Rj are the Carlson elliptic integrals of the 
second and third kind, available in the standard numerical 
libraries. 



S(p) = 



w(r)dr 



(B13) 



then the potential in the disc is 



V ( R >ti = yfi^ / ™M[(l+T)( P 2 +T)]*dT 

o, o — - 



(u+p,+T) 2 du 



1 / 2 [{U+P 2 +t)( U +1 + t)]2+2 



when o/O, and 
V(R,p) =AlnR + B(p,), 
with 



(B14) 



(B15) 



.1 = l(7E„ / ^(r)i? F (r + l,r + p 2 ,0)dr, 
ln(u+ ]1-\-t)w(t)AtAu 





oo oo 



B{p) = GE 



i 1 / 2 (u+p 2 + r) 1 / 2 (u+l + T)i 



(B16) 
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